#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include "aceMesh_runtime_c.h"
#include "cholesky_aceMesh.h"
#include "slave.h"
__thread_local_fix int i,j,k,s_starti,s_endi,s_startj,s_endj,s_startk,s_endk,s_NB;
__thread_local volatile unsigned long get_reply,put_reply;
// extern int NB;
//_thread_local_fix double A_slave[2][BLKX+2][BLKY+2][BLKZ+2];
__thread_local double s_T[_NB][_NB];
__thread_local double s_B[_NB][_NB];
__thread_local double s_C[_NB][_NB];
extern double A[_DIM][_DIM][_NB*_NB];

#define real double
#define integer int
int inline dgemm2(char *transa, char *transb, integer *m, integer *n, integer *k, real *alpha,
                        real *a, integer *lda, real *b, integer *ldb, real *beta, real *c__, integer *ldc)
{
    integer a_dim1, a_offset, b_dim1, b_offset, c_dim1, c_offset, i__1, i__2,i__3;
    real temp;
    integer i__, j, l;

    a_dim1 = *lda;
    a_offset = 1 + a_dim1 * 1;
    a -= a_offset;
    b_dim1 = *ldb;
    b_offset = 1 + b_dim1 * 1;
    b -= b_offset;
    c_dim1 = *ldc;
    c_offset = 1 + c_dim1 * 1;
    c__ -= c_offset;

    i__1 = *n;
    for (j = 1; j <= i__1; ++j)
    {
        i__2 = *k;
        for (l = 1; l <= i__2; ++l)
        {
            temp = *alpha * b[j + l * b_dim1];
            i__3 = *m;
            for (i__ = 1; i__ <= i__3; ++i__)
                c__[i__ + j * c_dim1] += temp * a[i__ + l * a_dim1];
        }
    }
    return 0;
}

int inline dtrsm2(char *side, char *uplo, char *transa, char *diag, integer *m,
                        integer *n, real *alpha, real *a, integer *lda, real *b, integer *ldb)
{
    integer a_dim1, a_offset, b_dim1, b_offset, i__1, i__2, i__3;
    real temp;
    integer i__, j, k;

    a_dim1 = *lda;
    a_offset = 1 + a_dim1 * 1;
    a -= a_offset;
    b_dim1 = *ldb;
    b_offset = 1 + b_dim1 * 1;
    b -= b_offset;

    i__1 = *n;
    for (k = 1; k <= i__1; ++k)
    {
        temp = 1.f / a[k + k * a_dim1];
        i__2 = *m;
        for (i__ = 1; i__ <= i__2; ++i__)
            b[i__ + k * b_dim1] = temp * b[i__ + k * b_dim1];
        i__2 = *n;
        for (j = k + 1; j <= i__2; ++j)
        {
            temp = a[j + k * a_dim1];
            i__3 = *m;
            for (i__ = 1; i__ <= i__3; ++i__)
                b[i__ + j * b_dim1] -= temp * b[i__ + k * b_dim1];
        }
    }
    return 0;
}

int inline dsyrk2(char *uplo, char *trans, integer *n, integer *k, real *alpha,
                        real *a, integer *lda, real *beta, real *c__, integer *ldc)
{
    integer a_dim1, a_offset, c_dim1, c_offset, i__1, i__2, i__3;
    real temp;
    integer i__, j, l;

    a_dim1 = *lda;
    a_offset = 1 + a_dim1 * 1;
    a -= a_offset;
    c_dim1 = *ldc;
    c_offset = 1 + c_dim1 * 1;
    c__ -= c_offset;

    i__1 = *n;
    for (j = 1; j <= i__1; ++j)
    {
        i__2 = *k;
        for (l = 1; l <= i__2; ++l)
        {
            temp = *alpha * a[j + l * a_dim1];
            i__3 = *n;
            for (i__ = j; i__ <= i__3; ++i__)
                c__[i__ + j * c_dim1] += temp * a[i__ + l * a_dim1];
        }
    }
    return 0;
}
#undef real
#undef integer

// -------------------------------------------------------------------------

void dpo(double *T, int N)
{
    int i, j, k;
    for(k = 0; k < N; k++)
    {
        double sum = 0;
        for(i = 0; i < k; i++)
            sum += T[i*N+k] * T[i*N+k];
        sum = T[k*N+k] - sum;
        T[k*N+k] = sqrt(sum > 0 ? sum : 0);
        for(i = k + 1; i < N; i++)
        {
            sum = 0;
            for(j = 0; j < k; j++)
                sum += T[j*N+i] * T[j*N+k];
            T[k*N+i] = (T[i*N+k] - sum) / T[k*N+k];
        }
    }

}


#define real double
#define integer int
int myspotf2_(char *uplo, integer *n, real *a, integer *lda, integer *info)
{
    static integer c__1 = 1;
    static real c_b10 = -1.f;
    static real c_b12 = 1.f;
    integer a_dim1, a_offset, i__1, i__2, i__3;
    real r__1;

    static integer j;
    static real ajj;
    extern double mysdot_(integer *, real *, integer *, real *, integer *);
    extern int mysscal_(integer *, real *, real *, integer *),
           mysgemv_(char *, integer *, integer *, real *, real *, integer *,
                    real *, integer *, real *, real *, integer *);
    a_dim1 = *lda;
    a_offset = 1 + a_dim1;
    a -= a_offset;
    i__1 = *n;
    for (j = 1; j <= i__1; ++j)
    {
        i__2 = j - 1;
        ajj = a[j + j * a_dim1] - mysdot_(&i__2, &a[j + a_dim1], lda, &a[j + a_dim1], lda);
        if (ajj <= 0.f)
        {
            a[j + j * a_dim1] = ajj;
            *info = j;
            return 0;
        }
        ajj = sqrt(ajj);
        a[j + j * a_dim1] = ajj;
        if (j < *n)
        {
            i__2 = *n - j;
            i__3 = j - 1;
            mysgemv_("No transpose", &i__2, &i__3, &c_b10, &a[j + 1 + a_dim1], lda, &a[j + a_dim1],
                     lda, &c_b12, &a[j + 1 + j * a_dim1], &c__1);
            i__2 = *n - j;
            r__1 = 1.f / ajj;
            mysscal_(&i__2, &r__1, &a[j + 1 + j * a_dim1], &c__1);
        }
    }
    return 0;
}

double mysdot_(integer *n, real *sx, integer *incx, real *sy, integer *incy)
{
    integer i__1;
    real ret_val;
    static integer i__, m, ix, iy, mp1;
    static real stemp;
    --sy;
    --sx;
    stemp = 0.f;
    ret_val = 0.f;
    if (*n <= 0)
    {
        return ret_val;
    }
    if (*incx == 1 && *incy == 1)
    {
        goto L20;
    }

    ix = 1;
    iy = 1;

    i__1 = *n;
    for (i__ = 1; i__ <= i__1; ++i__)
    {
        stemp += sx[ix] * sy[iy];
        ix += *incx;
        iy += *incy;
    }
    ret_val = stemp;
    return ret_val;

L20:
    m = *n % 5;
    if (m == 0)
    {
        goto L40;
    }
    i__1 = m;
    for (i__ = 1; i__ <= i__1; ++i__)
    {
        stemp += sx[i__] * sy[i__];
    }
    if (*n < 5)
    {
        ret_val = stemp;
        return ret_val;
    }

L40:
    mp1 = m + 1;
    i__1 = *n;
    for (i__ = mp1; i__ <= i__1; i__ += 5)
    {
        stemp = stemp + sx[i__] * sy[i__] + sx[i__ + 1] * sy[i__ + 1] + sx[i__ + 2] * sy[i__ + 2]
                + sx[i__ + 3] * sy[i__ + 3] + sx[i__ + 4] * sy[i__ + 4];
    }
}
int mysscal_(integer *n, real *sa, real *sx, integer *incx)
{
    integer i__1, i__2;
    static integer i__, m, mp1, nincx;
    --sx;
    if (*n <= 0 || *incx <= 0)
    {
        return 0;
    }
    if (*incx == 1)
    {
        goto L20;
    }
    nincx = *n **incx;
    i__1 = nincx;
    i__2 = *incx;
    for (i__ = 1; i__2 < 0 ? i__ >= i__1 : i__ <= i__1; i__ += i__2)
    {
        sx[i__] = *sa * sx[i__];
    }
    return 0;

L20:
    m = *n % 5;
    if (m == 0)
    {
        goto L40;
    }
    i__2 = m;
    for (i__ = 1; i__ <= i__2; ++i__)
    {
        sx[i__] = *sa * sx[i__];
    }
    if (*n < 5)
    {
        return 0;
    }

L40:
    mp1 = m + 1;
    i__2 = *n;
    for (i__ = mp1; i__ <= i__2; i__ += 5)
    {
        sx[i__] = *sa * sx[i__];
        sx[i__ + 1] = *sa * sx[i__ + 1];
        sx[i__ + 2] = *sa * sx[i__ + 2];
        sx[i__ + 3] = *sa * sx[i__ + 3];
        sx[i__ + 4] = *sa * sx[i__ + 4];
    }
    return 0;
}

int mysgemv_(char *trans, integer *m, integer *n, real *alpha, real *a,
             integer *lda, real *x, integer *incx, real *beta, real *y, integer *incy)
{
    integer a_dim1, a_offset, i__1, i__2;
    static integer i__, j, jx;
    static real temp;
    a_dim1 = *lda;
    a_offset = 1 + a_dim1;
    a -= a_offset;
    --x;
    --y;
    jx = 1;

    i__1 = *n;
    for (j = 1; j <= i__1; ++j)
    {
        if (x[jx] != 0.f)
        {
            temp = *alpha * x[jx];
            i__2 = *m;
            for (i__ = 1; i__ <= i__2; ++i__)
            {
                y[i__] += temp * a[i__ + j * a_dim1];
            }
        }
        jx += *incx;
    }
    return 0;
}
#undef real
#undef integer


void dpotrf_tile(double *T, int NB)
{
    int info;
    myspotf2_("L", &NB, T, &NB, &info);
}

void dgemm_tile(double *T, double *B, double *C, int NB)
{
    double done=1.0, dmone=-1.0;
    char TR='T', NT='N';
    dgemm2(&NT, &TR, &NB, &NB, &NB, &dmone, T, &NB, B, &NB, &done, C, &NB);
}

void dtrsm_tile(double *T, double *B, int NB)
{
    double done=1.0;
    char LO='L', TR='T', NU='N', RI='R';
    dtrsm2(&RI, &LO, &TR, &NU, &NB, &NB, &done, T, &NB, B, &NB);
}

void dsyrk_tile( double *T, double *B, int NB)
{
    double done=1.0, dmone=-1.0;
    char LO='L', NT='N';
    dsyrk2(&LO, &NT, &NB, &NB, &dmone, T, &NB, &done, B, &NB);
}
#ifdef MMEM
void slave_stencil_dgemm(targs1* args)
{
    dgemm_tile(args->T, args->B, args->C, args->NB);
}
 
void slave_stencil_dsyrk(targs2* args)
{
    dsyrk_tile(args->T, args->B, args->NB);
}

void slave_stencil_dpotrf(targs3* args)
{
    dpotrf_tile(args->T, args->NB);
}
                                                                                                                               
void slave_stencil_dtrsm(targs4* args)
{
    dtrsm_tile(args->T, args->B, args->NB);
}

#else
void slave_stencil_dgemm(targs1* args)
{
    s_starti=args->starti;
    s_startj=args->startj;
    s_startk=args->startk;
    get_reply=0;
    s_NB=args->NB;
    athread_get(PE_MODE,&A[s_startk][s_starti][0]
               , &s_T[0][0]
               , s_NB*s_NB*8, (void*)&get_reply
               , 0,0,0);
    athread_get(PE_MODE,&A[s_startk][s_startj][0]
               ,&s_B[0][0]
               ,s_NB*s_NB*8,(void*)&get_reply
               ,0,0,0);
    athread_get(PE_MODE,&A[s_startj][s_starti][0]
               ,&s_C[0][0]
               ,s_NB*s_NB*8,(void*)&get_reply
               ,0,0,0);           
    while(get_reply!=3);

    dgemm_tile((double*)&s_T, (double*)&s_B, (double*)&s_C, args->NB);

    put_reply=0;
   
    athread_put(PE_MODE,&s_C[0][0]
               ,&A[s_startj][s_starti][0]
               ,s_NB*s_NB*8,(void*)&put_reply
               ,0,0);           
    while(put_reply!=1);

}

void slave_stencil_dsyrk(targs2* args)
{
    s_starti=args->starti;
    s_startj=args->startj;
    get_reply=0;
    s_NB=args->NB;
    athread_get(PE_MODE,&A[s_starti][s_startj][0]
               , &s_T[0][0]
               , s_NB*s_NB*8, (void*)&get_reply
               , 0,0,0);
    athread_get(PE_MODE,&A[s_startj][s_startj][0]
               ,&s_B[0][0]
               ,s_NB*s_NB*8,(void*)&get_reply
               ,0,0,0);
    while(get_reply!=2);           
    dsyrk_tile((double*)&s_T,(double*)&s_B, args->NB);

    put_reply=0;
   

    athread_put(PE_MODE,&s_B[0][0]
               ,&A[s_startj][s_startj][0]
               ,s_NB*s_NB*8,(void*)&put_reply
               ,0,0);           
    while(put_reply!=1);

}

void slave_stencil_dpotrf(targs3* args)
{
    s_startj=args->startj;
    get_reply=0;
    s_NB=args->NB;
    athread_get(PE_MODE,&A[s_startj][s_startj][0]
               , &s_T[0][0]
               , s_NB*s_NB*8, (void*)&get_reply
               , 0,0,0);
    while(get_reply!=1);         
    dpotrf_tile((double*)&s_T, args->NB);
    put_reply=0;
    athread_put(PE_MODE,&s_T[0][0]
               ,&A[s_startj][s_startj][0]
               ,s_NB*s_NB*8,(void*)&put_reply
               ,0,0); 
    while(put_reply!=1);  
}

void slave_stencil_dtrsm(targs4* args)
{
    s_starti=args->starti;
    s_startj=args->startj;
    get_reply=0;
    s_NB=args->NB;
    athread_get(PE_MODE,&A[s_startj][s_startj][0]
               , &s_T[0][0]
               , s_NB*s_NB*8, (void*)&get_reply
               , 0,0,0);
    athread_get(PE_MODE,&A[s_startj][s_starti][0]
               ,&s_B[0][0]
               ,s_NB*s_NB*8,(void*)&get_reply
               ,0,0,0);
    while(get_reply!=2);
    dtrsm_tile((double*)&s_T, (double*)&s_B, args->NB);
    put_reply=0;
    athread_put(PE_MODE,&s_B[0][0]
               ,&A[s_startj][s_starti][0]
               ,s_NB*s_NB*8,(void*)&put_reply
               ,0,0); 
    while(put_reply!=1);  
}
#endif